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By means of extensive equilibrium molecular dynamics simulations we have investigated, the behav- 
ior of the interfacial tension 7 of two immiscible symmetrical Lennard-Jones fluids. This quantity 

k T 

is studied as function of reduced temperature T* — -B — in the range 0.6 < T* < 3.0. We find 
that, unlike the monotonic decay obtained for the liquid-vapor interfacial tension, for the liquid- 
liquid interface, j(T) has a maximum at a specific temperature. We also investigate the effect that 
surfactant-like particles has on the thermodynamic as well as the structural properties of the liquid- 
liquid interface. It is found that 7 decays monotonically as the concentration of the surfactant-like 
particles increases. 



I. INTRODUCTION 

The investigation of interfacial tension behavior in 
complex fluids is of relevance from tlie theoretical as well 
as from the practical point of viewlil. In particular, the 
nature of the liquid- vapot. interface has been extensively 
studiedrboth analyticallycl and numerically, using Monte 
CarlotTEl and Molecular DynamicsQtl algorithms. On the 
other hand, the properties of the liquid-liquid interface, 
which are of importance in different biological, systemsu 
as well as in various technological applications^, have re- 
ceived significantly less attention. One reason may be the 
complexity of the topology of the phase diagram for liquid 
mixtures. It is known that the intermolecular properties 
of a single component fluid have no important effect on 
the topology of its thermodynamic phase diagram lead- 
ing to the concept of corresponding states. By contrast, 
for liquid mixtures the interplay between the molecular 
properties of the components as well as the type of the 
intermolecular potential between different species leads 
to a highly noittrivial topology of the thermodynamic 
phase diagramtilO Recently, the properties of liquid- 
liquid planar interfaces have begun to be studied by 
means of molecular dynamics simulations using Lennard- 
Jones (LJ) interactionsfcaiia, by Monte Carlo simulations 
of benzene-water interfaced and iiexanol-waterEZI as well 
as by density functional theoryta. In addition, there 
have also been molecular dynamics simulations for the 



interface of water-alcohols using more complicated inter- 
molecular potentials-.that include the relevant molecular 
degrees of freedom.Ej~E3. The density and pressure pro- 
files of a symmetrical binary mixture of LJ fluids at the 
reduced temperature T* — 0.827 was studiecO in a re- 
gion of the phase diagram where the system is separated 
into two phases. It was found that as the system be- 
comes more miscible there is more diffusion of particles 
from one phase to the other and, as a consequence, the 
interfacial tension decreases. Some dynamic properties of 
this system were investigatedcJ at two reduced densities, 
p* =0.8 and 1.37 for several temperatures from which 
a crude sketch of the phase diagram was suggested. A 
third type of particle representing a simple model of am- 
phiphilic molecule was introduced finding that for low 
concentrations (of the order of 6%) there is no important 
effect on the binary fluid interface. In recent paperai3'E-3 
the structural properties of the above mentioned binary 
mixture was investigated at low temperatures and high 
pressures finding a stable oscillatory behavior in the den- 
sity profiles close to the interface. These oscillations de- 
crease as the temperature increases. It is also suggested 
that as pressure grows the interfacial tension increases. 
A recent-density functional investigation of a binary fluid 
mixtures shows that the interfacial tension as a function 
of temperature can have a maximum. 

With the purpose of getting a more complete under- 
standing of the thermodynamic and structural properties 
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of this liquid- liquid interface we have performed extensive 
equilibrium molecular dynamics simulations. The main 
results found in this work are: (i) the non-monotonic 
behavior of the interfacial tension as a function of tem- 
perature and (ii) the monotonic decay of this quantity as 
the concentration of a third type of particle, amphiphilic- 
like, increases. The layout of this paper is as follows: In 
section II we introduce the model and details of the sim- 
ulations. In section III we present the results for the 
thermodynamic as well as structural properties of the 
binary and ternary systems. Finally, in section IV the 
conclusions of this investigation are presented. 



II. MODEL AND SIMULATIONS. 

We consider a symmetrical binary mixture of partially 
miscible Lennard-Jones fluids. The intermolecular inter- 
actions F AA , F BB and F AB between particles A-A, B-B 
and A-B are described by modified LJ potentials that 
yield the forces, 
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where e and a are the same for all interactions and 
the parameter a XY controls the miscibility of the two 
fluids. For a partially miscible binary fluid we choose 
a AA = a BB = 1 and < a AB < 1. The interaction 
between particles of the same type is energetically more 
favourable, stronger binding, than the interaction A-B 
between different particles. For this reason one would 
expect a separation of the species at lower temperture 
when the entropy looses versus potential energy. With 
the aim of understanding the process by which a "sur- 
factant" weakens the interfacial tension and leads to the 
mixed state, a third species C is introduced to try to emu- 
late a surfactant-like particle. In this case we consider the 
same parameters in Eq. [l] and a cc = a AC = a BC = 1. 

We have carried out extensive equilibrium molecular 
dynamics (MD) simulations using the (N,V,T) ensemble 
for the particular miscibility parameter value a AB = 0.5. 
The range of the inter-molecular potential was set equal 
to three times the particle diameter a unless otherwise 
stated. In most of the simulations a total of N=1728 
particles were used. Nonetheless, to check for finite size 
effects we also carried out simulations with N=2592 par- 
ticles. These particles are placed in a paralelepiped of 
volume L x x L y x L Zl with L x = L y and L z = 2L Xl 
applying periodic boundary conditions in the x, y and z 
directions. The particles were initially placed on the sites 
of an FCC lattice forming a perfect planar interface, that 
is, all particles of type A are on the left side of the box 
while those of type B are on the opposite side. In this 
way one can obtain a minimum of two interfaces due to 
periodic boundary conditions. If one starts with a sta- 
tistical mixture usually more than two demixed regions 
develop giving rise to more than two interfaces. 



The three component system initial configuration was 
chosen from an equilibrated separated binary system 
picking at random N c /2 particles from type A and N c /2 
particles from type B and replacing them by particles 
of type C. This way of putting the third species in the 
system makes the total density to remain constant. It is 
customary to carry out the simulations using the follow- 
ing reduced units for the distance r* = — , particle linear 
momentum p* = , p and time t* = —\/e/m. In these 
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definitions m is the mass of each particle, which is taken 
to be the same for all particles, a is the particle diam- 
eter and e is the depth of the LJ potential. Similarly, 
one can define reduced thermodynamic quantities as fol- 
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lows: T* = — ^— that represents the reduced temperature 
with fc B =Boltzmann's constant and p* = pa 3 for the re- 
duced density, with p — N/V. The equations of motion 
were integrated using a fourth order predictor-corrector 
algorithm with an integration step-size of At* < 0.005, 
which in standard units is of the order of I0 -5 nanosec- 
onds in the scale of argon. The particles initial velocities 
were assigned from a Boltzmann distribution. The equi- 
libration times for most of the simulations were of the 
order of 10 time steps. Thermodynamic quantities were 
measured every 50 time-step iterations up to a total of 
5xl0 5 to I0 6 measurements from which averages were 
evaluated. This amounts to a simulation time between 
5-10 nanoseconds in the scale of argon. At the start the 
reduced homogeneous density in the simulation cell was 
set equal to p* = 0.844, which is close to density of the 
triple point of argon. Due to the inhomogeneity which 
develops in the system around the interface, the densi- 
ties of the bulk phases later are slightly higher than this 
starting density. The bulk densities are evaluated when 
the system has been equilibrated and they depend on the 
conditions of the system and on T* . 



III. RESULT AND DISCUSSION 

A. Binary mixture 

Since our interest is to study the structural and ther- 
modynamic properties of the interface, we have carried 
out a set of simulations for a sequence of temperatures 
below the critical demixing temperature for two system 
sizes, namely, N = 1728 and N = 2592 particles. We 
have considered these values of N to find out about pos- 
sible finite size effects. All the quantities studied show 
qualitatively the same tendency for the two system sizes. 
The interfacial behavior is investigated by calculating the 
density profiles, the pressure, and the interfacial tension 
at different temperatures. The relevant parameters of 
the investigated systems are summarized in Tables [j| and 
U where the bulk-values of the density and the pressure 
are specified. To emphasize the separated nature of the 
system, the values of the reduced total density in the 
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Sys. 


T* 


p bu.k 


„bulk A 

Pa 


bulk A 
Pb 


p* 


7* 


1 


0.6 


0.888 


0.888 


0.000 


0.587 


2.068±0.012 


2 


0.827 


0.876 


0.876 


0.000 


1.964 


2.143±0.012 


3 


1.1 


0.870 


0.870 


0.000 


3.510 


2.186±0.015 


4 


1.4 


0.866 


0.866 


0.000 


5.081 


2.164±0.018 


5 


1.7 


0.863 


0.863 


0.000 


6.555 


2.092±0.021 


6 


2.1 


0.858 


0.858 


0.000 


8.404 


1.941±0.033 


7 


2.6 


0.860 


0.849 


0.011 


10.651 


1.689±0.039 


8 


3.0 


0.859 


0.842 


0.017 


12.368 


1.539±0.044 



TABLE I. Results for the reduced bulk densities and re- 
duced normal pressure (P* = — P n ) for the binary mixture 
with N = 1728 particles. 



Sys. T* p bulk p bulkA p bulkA P* -f_ 

1 0.6 0.882 (0.7) 0.882 0.000 0.458™(22) 1.986±0.020 (4) 

2 0.827 0.873 (0.3) 0.873 0.000 1.876 (4.5) 2.103±0.020 (1.9) 

3 1.1 0.869 (0.1) 0.869 0.000 3.456 (1.5) 2.180±0.030 (0.3) 

4 1.4 0.865 (0.1) 0.865 0.000 5.068 (0.3) 2.190±0.050 (1.2) 
6 2.1 0.859 (0.1) 0.859 0.000 8.439 (0.4) 1.909±0.040 (1.7) 
8 3.0 0.860 (0.1) 0.842 0.018 12.472(0.8) 1.511±0.060 (1.8) 



TABLE II. Results for the reduced bulk densities and re- 
duced normal pressure for the binary mixture with N = 2592 
particles. The numbers in parenthesis are the percent differ- 
ences calculated with respect to the corresponding values in 
table | 



bulk A-rich phase p bulkA ) the reduced density of parti- 
cles A, p h 1 ulkA , and particles B, / o b J ulkA , are given. Due 
to the symmetry of the interactions, the B-rich phase is 
symmetric to the A-rich phase. 

At sufficiently low temperatures the reduced density 
of B particles in the phase of A is very small indicat- 
ing that the system is fully separated. As temperature 
increases the total density of the bulk region decreases 
slightly because the diffusion of particles through the in- 
terface increases, making the inhomogeneity smaller and 
driving the system towards a mixed state. This can be 
seen in Figures [j] and || where the density profiles along 
the z direction (longer side of the box) are shown for the 
N = 1728 and 2592 systems, respectively. The contin- 
uous line corresponds to a temperature of T* — 0.827, 
which can be considered low, (systems 2 in Tables [| and 
|J) while the dashed line corresponds to a higher temper- 
ature T* = 1.4 (systems 4 in Tables | and 0). 

When looking at the density profiles for the two system 
sizes studied, Figures | and |, we observe that at T* = 
0.827 these quantities exhibit an oscillatory structure at 
the interfaces. This reminds us about the behavior of 
the density profile in front of a hard wall, as discussed in 
referenceLj. This kind of structure is also found in the 
liquid- vapor system, but the oscillations are significantly 
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FIG. 1. Density profiles for the separated binary mixture 
with N=1728 particles and an average reduced total density 
p* = 0.844. 



stronger in the liquid-liquid interface. As expected, one 
also sees from these figures that the oscillations in the 
bulk region are less pronounced for the larger system. 

The values of the reduced normal pressure P* shown 
in column sixth of Tables Q and || follow an almost linear 
behavior as a function of temperature in the range stud- 
ied. The normal and transversal pressures profiles P n (z), 
Pt(z) were calculated using the definition of the Irving- 
Kir kwoO|d_piressure tensor, which for a planar interface is 
given byOQ 

p n{z) = p{z)k B T 
Pt(z) = p{z)k B T 

In these expressions the first term corresponds to the 
ideal gas contribution while the second comes from the in- 
termolecular forces. For a dense system and long ranged 
intermolecular potentials the latter term yields the larger 
contribution. On the left side of Figure |^ we show the 
normal and tangential components of the pressure tensor 
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FIG. 2. Density profiles for the separated binary mixture 
with N=2592 particles and an average reduced total density 
p* = 0.844. 
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FIG. 3. Pressure profiles for the system with TV = 1728 par- 
ticles. Left: Normal and transversal reduced pressure profiles 
for the separated binary mixture. Right: Ideal gas contribu- 
tion (Pineal — p*T*) as well as contributions from the interac- 
tion part to the transversal (-P t c ) and normal (Pn) components 
of the pressure tensor. 

as functions of z for T* = 0.827 and T* = 1.4. The first 
feature to notice is that P„(z) remains constant for both 
temperatures as one would expect for a system in ther- 
modynamic equilibrium. The pressures P* range from 
0.25 kbar to 5.18 kbar in the scale of argon for the sys- 
tems of Tables | and E3. On the right side of the same 
figure the ideal as well as the configurational parts of 
P*(z) and Pt(z) are plotted. At the lower temperature 
T* = 0.827 the pressure profiles clearly show the oscil- 
lations related to the oscillatory behavior of the density 
profiles. One can learn that the configurational contribu- 
tion is much larger than the ideal gas part since we are 
dealing with a dense system. These profiles are needed 
to evaluate the inlarfacial tension by means of the me- 
chanical definitioro. 

/•bulks 

7= / [P n (z)-P t {z)]dz. (3) 

J bulkA 

Also, a straightforward evaluation of 7 can be done using 
the Kirkwood-Buff formulaEffl 

" 4z(g( ] "' [r )■ (4) 

In this latter equation there appears an additional fac- 
tor of i which comes from having two interfaces in the 
system. From the calculational point of view it is better 
to use the Kirkwood-Buff formula rather than the me- 
chanical definition since fluctuations in P n (z) and Pt(z) 
may introduce important inaccuracies in the evaluation 
of 7 according to Eq. |3j. As a matter of fact, we calcu- 
lated the interfacial tension using both expressions and 
found consistency in the values within the statistics of 
the simulations. The range of temperatures, in reduced 
units, in which 7 has been studied is 0.6 <T*< 3.0, that 
in the scale of argon corresponds to 72 K < T < 359.5 K. 

Unlike for the liquid-vapor interface we find a non- 
monotonic behavior of j(T) for the liquid- liquid inter- 
face. In Figure we show the behavior of 7* as a function 



FIG. 4. Behavior of the reduced interfacial tension 
(7* = —7) as a function of reduced temperature for the sys- 
tems with TV = 1728 (•) and TV = 2592 (o) particles. All lines 
are a guide to the eye. The explicit values of 7* together with 
other thermodynamic variables are given in tables §andg 

of T* for systems with TV = 1728 and TV = 2592 parti- 
cles. The main feature of this graph is the maximum of 
the interfacial tension at about T* w 1.1. Such a maxi- 
mum has been reported recently in the context of density 
functional theoryla. It can be understood physically as 
follows: A stronger mixing near the interface introduces 
weaker A-B bonds and raises the potential energy. This 
leads to an increase of 7, the free energy of the inter- 
face, until at high temperatures the entropy contribution 
drags it down. We have calculated j(T) for these two 
systems to chek for possible finite size effects. In doing 
so we have to assure that the bulk thermodynamic quan- 
tities are approximately the same for both systems. In 
particular, we monitored very closely the bulk density 
since small variations in this quantity lead to significant 
changes in the pressure of the system at low temper- 
atures. So, the important key is to have approximately 
the same pressure for both systems. This can be achieved 
by adjusting the side L z of the simulational box for the 
larger system, mantaining its cross section constant. In 
this way we eliminate possible variations of j(T) due to 
variations in the cross section area. We found that the 
optimal average value of LI was 30 for the temperature 
range studied. In Figure f| and tables | and [H] one also 
sees that at low temperatures, the values of 7*(T*) (be- 
tween the two system sizes) are at the most 4% different, 
while at higher temperatures (T* > 1.1) the values are 
much closer to each other. One also observes that small 
variations in the bulk density lead to important changes 
in the pressure at low temperatures. However, at higher 
temperatures the variations in pressure are smaller. This 
is the expected behavior for a single LJ fluid. Therefore, 
differences in j(T) are due to the variations in pressure, 
although the former are much smaller than the latter. 
This happens because the interfacial tension is obtained 
from the difference between the normal and tangential 
components of the pressure tensor. This behavior ex- 
plains why the differences and consistencies in the values 
of 7*(T*). 
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FIG. 5. Density profiles for the three component system for 
different values of TV C at T* = 1.4 and total average reduced 
density p* = 0.844. The total number of particles in the 
system is TV = 1728. Continuous line correspond to species 
A, the long dashed to species B, dashed line to species C, and 
dotted line to the total density profile. 



B. Ternary mixture 

In this section we investigate the role that a third 
species plays in the interface properties of a demixed 
binary system. A third species — C-particles — is in- 
troduced in the system and the behavior of the inter- 
facial tension is studied as a function of concentration. 
The interactions of the C-particles F cc = F CA = F CB 
are all assumed equal to the strong interactions F AA — 
F BB = F cc . When placed between the A and B, C par- 
ticles avoid the weak A-B bonds and lower the potential 
energy. It is found that as N c increases 7 decays mono- 
tonically. This is reminiscent of surfactant-like behavior 
in some ternary systems. In Figure [| we plot the density 
profiles of the ternary system as well as the reduced to- 
tal density for different concentrations of C particles at 
T* = 1.4. These profiles yield a clear evidence that C- 
particles like to be at the interface position rather than 
in the bulk phases. In this way the C-particles dimin- 
ish (screen) the energetically unfavourable interactions 
between particles A-B. In Table III a summary of the re- 
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FIG. 6. Reduced interfacial tension as a function of iV c for 
a system with TV = 1728 particles. The inset shows the excess 
T c of particles C at the interface as a function of TV C . 



TABLE III. Results for the reduced bulk densities for the 
ternary mixture. 



duced bulk densities of particles A, B and C in the A-rich 
phase for T* = 0.827 and T* = 1.4 is given. As TV C in- 
creases the density of B in the A-rich phase also increases. 
This means that the C particles help the B particles to 
diffuse into the A-phase because the B particles can be 
solvated by C and therefore avoid the weak A-B interac- 
tions. As TV C becomes sufficiently large this mechanism 
drives the system to a mixed state. Therefore, the C par- 
ticles act like a surfactant or emulgator. This mechanism 
may have potential technological implications when try- 
ing to design compatibilizers. In Figure ^ the reduced 
interfacial tension is plotted as a function of TV C for two 
different temperatures. With increasing bulk concentra- 
tions of C-particles we see a monotonic decay. This low 
concentration behavior is consistent with that found in 
referenced where amphiphile like-particles C are intro- 
duced. In such a model 7* as a function of TV C shows a 
linear decay in the small amphiphile concentration range. 

In this same figure we also observe that the curve 
7*(iV c ) at T* = 1.4 is slightly above the curve of 7*(TV C ) 
at T* — 0.827. This is an unexpected result since usually 
7* decreases when T* increases. Nonetheless, we should 
recall that these two points lie in the region where 7* has 
its maximum as a function of T* for the binary system. 
This fact might also be responsible for the shoulder in the 
curve for 7* at T* = 1.4 in the region 200 < TV C < 300. 
We would like to emphasize that this structure is outside 
of the statistics of the simulations since the error bars 
are smaller than the difference between the points. We 
have also evaluated the excess r c of C-particles given by 
Equation || and plot this quantity as a function of TV C in 
the inset of Figure ^|. 
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FIG. 7. The reduced normal pressure P* = Pa 3 /e as a 
function of JV C at two different temperatures for a system 
with N = 1728 particles. 
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We find that for all values of N c , T C (T* = 0.827) > 
T C (T* = 1.4) which means that there are more parti- 
cles of type C located at the interface at the lower tem- 
perature and therefore they screen more effectively the 
energetically unfavourable interactions between A and B 
particles, diminishing even more the interfacial tension. 
In the range 200 < N c < 400 the excess of C parti- 
cles at the interface at T* — 1.4 increases slower than at 
T* = 0.827 giving rise to a shoulder in 7* as shown in 
figure [| 

Since in the above results not the pressure but the over- 
all density p* — 0.844 was kept constant, in Figure |t] we 
plot the behavior of P* as a function of N c for the two 
temperatures indicated above. One notes that the hydro- 
static pressure decreases slowly as N c increases. This is 
due to the fact that N c /2 particles from fluid A and the 
same amount from fluid B were switched to particles C, 
replacing the weak A-B attraction by the stronger C-A 
and C-B attraction. This increases the cohesion and re- 
duces the pressure. One could be tempted to obtain the 
upper curve from the lower by simply adding the ideal 
gas term p*AT*, however, this is not the case because, 
as noted above with regard to Figure the contribution 
of the configurational part to the pressure is the more 
relevant one. 



IV. CONCLUSIONS 

Molecular dynamics simulations for a mixture of two 
kinds of molecules A and B lead to a phase separation 
with an A-rich and B-rich phase and an interface between 
the phases. This happens if the attractive potential be- 
tween A-B is weaker than that between A-A and B-B and 
if density and temperature are in the two phase region 
of demixing. The interface between the demixed phases 
has also been investigated. We have calculated density 
profiles at different temperatures, the components of the 



pressure tensor across the interface and have evaluated 
the interfacial tension 7 for a series of temperatures. The 
important result of the first part of this investigation is 
the maximum in j(T) (Figure ^). We explain this result 
as a consequence of a strong mixing near the interface 
that yield weaker A-B bonds thus raising the potential 
energy, that in turn leads to an increase in j{T), until at 
higher temperatures the interfacial free energy decreases 
due to an increase of the entropy term. 

Then we have added a third kind of particle, 
"surfactant- like" , and have studied the interfacial tension 
as a function of its concentration. We found that these 
particles reduce the interfacial tension 7 as its concen- 
tration increases. This happens because these particles 
screen the unfavorable A-B interactions and lead to a 
mixing state of the system. 
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